Pylops - acoustic wavefield separation

Author: M.Ravasi

This notebook is inspired by the paper of J. van der Neut and F. J. Herrmann, Up/Down Wavefield Decomposition by Sparse Inversion, 74th EAGE Conference and Exhibition, where multi-component seismic data ($p(t,x)$ and $v_z(t,x)$) are inverted by inverting the underlying equations that explain those wavefields as superposition of the up- and downgoing components of the pressure wavefield ($p^-(t,x)$ and $p^+(t,x)$) :

$$ \begin{bmatrix} \mathbf{p} \\ \mathbf{v_z} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ \frac{k_z}{\omega \rho} & - \frac{k_z}{\omega \rho} \\ \end{bmatrix} \begin{bmatrix} \mathbf{p^+} \\ \mathbf{p^-} \end{bmatrix} = \mathbf{W} \mathbf{p_{dec}}$$

as in Wapenaar, K., 1998, Reciprocity properties of one-way propagators: Geophysics, Vol. 63, 1795-1798, where all wavefields are expressed here in the $(f, k_x)$ domain in 2d (or $(f, k_x, k_y)$ in 3d), $\omega=2\pi f$ is the angular frequency, $k_x$ and $k_y$ are the horizontal wavenumbers, and $k_z=\sqrt{k^2-k_x^2}$ is the vertical wavenumber. Morevoer, $\frac{k_z}{\omega \rho}$ is called obliquity factor.

To start we will consider the analytic solution of the above equation as in practice we record $p$ and $v_z$.

$$ \begin{bmatrix} \mathbf{p^+} \\ \mathbf{p^-} \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 1 & \frac{\omega \rho}{k_z} \\ 1 & - \frac{\omega \rho}{k_z} \\ \end{bmatrix} \begin{bmatrix} \mathbf{p} \\ \mathbf{v_z} \end{bmatrix}$$

The advantage of obtaining $p^-$ and $p^+$ by inversion is that we can stabilize our procedure and obtained improved wavefields in presence of noisy data as shown in van der Neut and Herrman (2015). Moreover as the equations are in the wavenumber domain, having densely, regularly sampled x and y axes (receivers in IL and XL) direction is a pre-requisite for performing the FFT along those axes (while it is possible to perform FFT with irregularly sampled axis, this is costly, and sparse sampling leads to aliasing which does not allow appyling the obliquity factor to $V_z$).

A possible way to circumvent this condition (when not achieved in acquisition) is to write the problem in the following way:

$$ \mathbf{d} = \frac{1}{2} \mathbf{R} \mathbf{F}^H \mathbf{W} \mathbf{F} \mathbf{p_{dec}} $$

where $\mathbf{F}$ is a 2d or 3d FFT, $\mathbf{R}$ is a restriction (or even better a bilinear interpolation) operator, and the model and data are now in $(t,x)$ domain. The actual wavefield combination via $\mathbf{W}$ is done on regularly sampled wavefields and then sampling in performed in the original time-space domain to bring them to the actual acqisition grid. When solved by inversion this problem allows us to retreive up- and down- decomposed field in an ideal, regualrly and densely sampled geometry. It is expected that to be able to solve this problem a sparsity-promoting condition should be defined for the model to retrieved model such as

$$J = |\mathbf{x_{dec}}|_1 s.t. \mathbf{d} = \frac{1}{2} \mathbf{R} \mathbf{F}^H \mathbf{W} \mathbf{F} \mathbf{S}^H \mathbf{x_{dec}} $$

where $\mathbf{x_{dec}}=\mathbf{S}^H \mathbf{x_{dec}}$ are the wavefields in a sparse domain (ideally curvelets or wave atoms). Note that if f-k can be the sparse domain then the problem becomes even simpler

$$J = |\mathbf{x_{dec}}|_1 s.t. \mathbf{d} = \frac{1}{2} \mathbf{R} \mathbf{F}^H \mathbf{W} \mathbf{x_{dec}} $$

This notebooks is organized as follow:

  1. 2D analytical wavefield separation using PyLops chained operators
  2. 2D wavefield separation by inversion
  3. 2D analytical wavefield separation with irregularly and sparsely sampled x and sparsity promoting inversion (TO BE DONE..)
  4. 3D analytical wavefield separation by inversion (TO BE DONE..)
  5. 3D analytical wavefield separation with irregularly and sparsely sampled x and y and sparsity promoting inversion (TO BE DONE..)
In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp

from scipy.sparse import csr_matrix, vstack
from scipy.signal import filtfilt
from scipy.linalg import lstsq, solve
from scipy.sparse.linalg import LinearOperator, cg, lsqr
from scipy import misc

from pylops.utils                      import dottest
from pylops.utils.wavelets             import *
from pylops.utils.seismicevents        import *
from pylops.basicoperators             import *
from pylops.signalprocessing import *
from pylops.waveeqprocessing.wavedecomposition import *


from pylops.optimization.leastsquares  import *
from pylops.optimization.sparsity  import *

Let's import input data from fiel

In [2]:
inputfile = '../../pylops/testdata/updown/input.npz'
vel_sep = 2400.0         # velocity at separation level
rho_sep = 1000.0         # density at separation level

inputdata = np.load(inputfile)

# Receivers
r = inputdata['r']
nr = r.shape[1]
dr = r[0, 1]-r[0, 0]

# Sources
s = inputdata['s']

# Model
rho = inputdata['rho']

# Axes
t = inputdata['t']
nt, dt = len(t), t[1]-t[0]
x, z = inputdata['x'], inputdata['z']

# Data
p = inputdata['p'].T
vz = -inputdata['vz'].T

# Normalize
vz /= p.max()
p /= p.max()

# Add noise
#p += np.random.normal(0, 1e4, p.shape)
#vz += np.random.normal(0, 2e-2, p.shape)

# First arrival
direct = np.sqrt(np.sum((s[:,np.newaxis]-r)**2, axis=0))/vel_sep

plt.figure(figsize=(9,5))
plt.imshow(rho, cmap='gray', extent = (x[0], x[-1], z[-1], z[0]))
plt.scatter(r[0, ::5],  r[1, ::5], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(s[0], s[1], marker='*', s=250, c='r', edgecolors='k')
plt.axis('tight')
plt.xlabel('x [m]'),plt.ylabel('y [m]'),plt.title('Model and Geometry')
plt.xlim(x[0], x[-1]);
plt.ylim(z[-1], z[0]);

fig, axs = plt.subplots(1, 2, figsize=(15,12))
im=axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$P$')
axs[0].axis('tight')
plt.colorbar(im, ax=axs[0])
im=axs[1].imshow(vz.T, cmap='gray', vmin=-0.1*np.abs(vz).max(), vmax=0.1*np.abs(vz).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$V_z$')
axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
plt.colorbar(im, ax=axs[1]);

Create frequency-wavenumber spectra (using PyLops FFT2D)

In [3]:
nfft=2**11
FFTop = FFT2D(dims=[nr, nt], 
              nffts=[nfft, nfft], sampling=[dr, dt])
dottest(FFTop, nfft*nfft, nt*nr, complexflag=2)

P = FFTop*p.flatten()
VZ = FFTop*vz.flatten()
P = P.reshape(nfft, nfft)
VZ = VZ.reshape(nfft, nfft)

1a. 2D Analytical wavefield separation

In [4]:
critical = 1.1
ntaper = 51

#obliquity factor
[Kx, F] = np.meshgrid(FFTop.f1, FFTop.f2, indexing='ij')
k=F/vel_sep
#Kz=np.sqrt(k**2-Kx**2)
Kz=np.sqrt((k**2-Kx**2).astype(np.complex))
Kz[np.isnan(Kz)] = 0

OBL=rho_sep*(np.abs(F)/Kz)
OBL[Kz==0]=0

# cut off and taper
mask = np.abs(Kx)<critical*np.abs(F)/vel_sep
OBL = OBL*mask
OBL = filtfilt(np.ones(ntaper)/float(ntaper), 1, OBL, axis=0)
OBL = filtfilt(np.ones(ntaper)/float(ntaper), 1, OBL, axis=1)

OBL1 = Kz /(rho_sep*np.abs(F))
OBL1[F==0] = 0
OBL1 = OBL1*mask
OBL1 = filtfilt(np.ones(ntaper)/float(ntaper), 1, OBL1, axis=0)
OBL1 = filtfilt(np.ones(ntaper)/float(ntaper), 1, OBL1, axis=1)

# scaled Vz
VZ_obl = OBL*VZ;
vz_obl = FFTop.H*VZ_obl.flatten()
vz_obl = np.real(vz_obl.reshape(nr, nt))

p = FFTop.H*P.flatten()
p = np.real(p.reshape(nr, nt))

# separation
pup=(p-vz_obl)/2;
pdown=(p+vz_obl)/2;


fig, axs = plt.subplots(1, 4, figsize=(19,5))
im=axs[0].imshow(np.fft.fftshift(np.abs(P[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
             extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],),
             vmin=0, vmax=np.abs(P).max())
axs[0].set_title(r'$P(f, k_x)$')
axs[0].axis('tight')
plt.colorbar(im, ax=axs[0])
im=axs[1].imshow(np.fft.fftshift(np.abs(Kz[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
             extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],))
axs[1].set_title(r'$K_z$')
axs[1].axis('tight')
plt.colorbar(im, ax=axs[1])
im=axs[2].imshow(np.fft.fftshift(np.abs(OBL[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],))
axs[2].set_title(r'$OBL$')
axs[2].axis('tight')
plt.colorbar(im, ax=axs[2])
im=axs[3].imshow(np.fft.fftshift(np.abs(VZ_obl[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],),
                 vmin=0, vmax=np.abs(P).max())
axs[3].set_title(r'$VZ_obl$')
axs[3].axis('tight')
plt.colorbar(im, ax=axs[3])

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(vz_obl.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$vzobl$')
axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(pup.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
#axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$pup$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(pdown.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$pdown$')
#axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr//2], 'r', lw=2, label=r'$p$')
plt.plot(t, vz_obl[nr//2], '--b', lw=2, label=r'$v_z$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.figure(figsize=(14, 3))
plt.plot(t, pup[nr//2], 'r', lw=2, label=r'$p^-$')
plt.plot(t, pdown[nr//2], '--b', lw=2, label=r'$p^+$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
Out[4]:
<matplotlib.legend.Legend at 0x1c1a173470>

1b. Analytical wavefield separation with PyLops operators

In [5]:
OBLop = Diagonal(OBL.flatten(), dtype='complex128')

Sop = 0.5*(BlockDiag([FFTop.H, FFTop.H])*\
    Block([[Identity(nfft*nfft, dtype='complex128'), OBLop],
           [Identity(nfft*nfft, dtype='complex128'), -OBLop]])*\
           BlockDiag([FFTop, FFTop]))

d = np.concatenate((p.flatten(), vz.flatten()))
dud_lop = np.real(Sop*d)

d = d.reshape(2*nr, nt)
dud_lop = dud_lop.reshape(2*nr, nt)

pdown_lop, pup_lop= dud_lop[:nr], dud_lop[nr:]

print(np.linalg.norm(pdown - pdown_lop))
print(np.linalg.norm(pup - pup_lop))
1.401961443511688e-15
1.089710079892621e-15
In [6]:
fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
#axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(vz_obl.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$vzobl$')
#axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(pup_lop.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
#axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$pup$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(pdown_lop.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$pdown$')
#axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr//2], 'r', lw=2, label=r'$p$')
plt.plot(t, vz_obl[nr//2], '--b', lw=2, label=r'$v_z$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.figure(figsize=(14, 3))
plt.plot(t, pup_lop[nr//2], 'r', lw=2, label=r'$p^-$')
plt.plot(t, pdown_lop[nr//2], '--b', lw=2, label=r'$p^+$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend();
In [7]:
PUP_lop = (FFTop*pup_lop.flatten()).reshape(nfft, nfft)
PDOWN_lop = (FFTop*pdown_lop.flatten()).reshape(nfft, nfft)

fig, axs = plt.subplots(1, 2, figsize=(10,5))
im=axs[0].imshow(np.fft.fftshift(np.abs(PDOWN_lop[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],),
                 vmin=0, vmax=np.abs(P).max())
axs[0].set_title(r'$P^+(f, k_x)$')
axs[0].axis('tight')
axs[0].set_ylim(80, 0);
im=axs[1].imshow(np.fft.fftshift(np.abs(PUP_lop[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],),
                 vmin=0, vmax=np.abs(P).max())
axs[1].set_title(r'$P^-(f, k_x)$')
axs[1].axis('tight');
axs[1].set_ylim(80, 0);

2. Wavefield separation by inversion

In [8]:
OBLop = Diagonal(OBL1.flatten(), dtype='complex128')

S1op_scaled = (BlockDiag([FFTop.H, (p.max()/vz.max())*FFTop.H])*\
        Block([[Identity(nfft*nfft, dtype='complex128'), Identity(nfft*nfft, dtype='complex128')],
               [OBLop, -OBLop]])*\
        BlockDiag([FFTop, FFTop]))

print(p.max()/vz.max())
d = np.concatenate((p.flatten(), (p.max()/vz.max())*vz.flatten()))
#dottest(S1op_scaled, 2*nt*nr, 2*nt*nr, tol=1e0)
2249818.971759595
In [9]:
dud_inv, istop, itn, r1norm, r2norm = \
    lsqr(S1op_scaled, d.flatten(), damp=1e-10, 
         iter_lim=20, show=2)[0:5]
dud_inv = np.real(dud_inv)
dud_inv = dud_inv.reshape(2*nr, nt)

pdown_inv, pup_inv= dud_inv[:nr], dud_inv[nr:]

print(np.linalg.norm(pdown - pdown_inv))
print(np.linalg.norm(pup - pup_inv))

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
#axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(vz_obl.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$vzobl$')
#axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(pup_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
#axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$pup$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(pdown_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$pdown$')
#axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr//2], 'r', lw=2, label=r'$p$')
plt.plot(t, vz_obl[nr//2], '--b', lw=2, label=r'$v_z$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.figure(figsize=(14, 3))
plt.plot(t, pup_inv[nr//2], 'r', lw=2, label=r'$p^-$')
plt.plot(t, pdown_inv[nr//2], '--b', lw=2, label=r'$p^+$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   434562 rows  and   434562 cols
damp = 1.00000000000000e-10   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       20
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   1.361e+01  1.361e+01    1.0e+00  9.8e-02
     1  3.72488e-06   3.221e+00  3.221e+00    2.4e-01  6.4e-01   1.4e+00  1.0e+00
     2 -4.90167e-06   1.643e+00  1.643e+00    1.2e-01  4.2e-01   1.7e+00  2.1e+00
     3 -4.90827e-06   9.905e-01  9.905e-01    7.3e-02  3.2e-01   2.0e+00  3.4e+00
     4 -2.74583e-06   6.425e-01  6.425e-01    4.7e-02  2.7e-01   2.2e+00  4.9e+00
     5 -7.91046e-06   4.318e-01  4.318e-01    3.2e-02  2.4e-01   2.4e+00  6.4e+00
     6 -2.99539e-06   2.958e-01  2.958e-01    2.2e-02  2.2e-01   2.6e+00  8.1e+00
     7 -6.02954e-06   2.087e-01  2.087e-01    1.5e-02  1.9e-01   2.8e+00  9.7e+00
     8 -5.05354e-06   1.549e-01  1.549e-01    1.1e-02  1.6e-01   2.9e+00  1.1e+01
     9 -6.40115e-06   1.227e-01  1.227e-01    9.0e-03  1.4e-01   3.1e+00  1.3e+01
    10 -5.48717e-06   1.034e-01  1.034e-01    7.6e-03  1.1e-01   3.2e+00  1.5e+01
    11 -6.22387e-06   9.155e-02  9.155e-02    6.7e-03  8.8e-02   3.4e+00  1.7e+01
    12 -5.83974e-06   8.401e-02  8.401e-02    6.2e-03  7.3e-02   3.5e+00  1.8e+01
    13 -6.40602e-06   7.881e-02  7.881e-02    5.8e-03  6.2e-02   3.6e+00  2.1e+01
    14 -6.49674e-06   7.484e-02  7.484e-02    5.5e-03  5.5e-02   3.7e+00  2.3e+01
    15 -6.16872e-06   7.153e-02  7.153e-02    5.3e-03  5.0e-02   3.8e+00  2.5e+01
    16 -6.73433e-06   6.866e-02  6.866e-02    5.0e-03  4.7e-02   4.0e+00  2.8e+01
    17 -6.16827e-06   6.617e-02  6.617e-02    4.9e-03  4.3e-02   4.1e+00  3.1e+01
    18 -6.84475e-06   6.406e-02  6.406e-02    4.7e-03  3.9e-02   4.2e+00  3.4e+01
    19 -6.63506e-06   6.232e-02  6.232e-02    4.6e-03  3.5e-02   4.3e+00  3.7e+01
    20 -6.64805e-06   6.087e-02  6.087e-02    4.5e-03  3.1e-02   4.4e+00  4.0e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 6.1e-02   anorm = 4.4e+00   arnorm = 8.4e-03
itn   =      20   r2norm = 6.1e-02   acond = 4.0e+01   xnorm  = 1.1e+01
 
0.877757730395036
0.8777577309633208
Out[9]:
<matplotlib.legend.Legend at 0x1c19097d30>
In [10]:
PUP_inv = (FFTop*pup_inv.flatten()).reshape(nfft, nfft)
PDOWN_inv = (FFTop*pdown_inv.flatten()).reshape(nfft, nfft)

fig, axs = plt.subplots(1, 2, figsize=(19,5))
im=axs[0].imshow(np.fft.fftshift(np.abs(PDOWN_inv[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],),
                 vmin=0, vmax=np.abs(P).max())
axs[0].set_title(r'$P^+(f, k_x)$')
axs[0].axis('tight')
axs[0].set_ylim(80, 0);
plt.colorbar(im, ax=axs[0])
im=axs[1].imshow(np.fft.fftshift(np.abs(PUP_inv[:, :nfft//2-1]),axes=0).T, cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1], FFTop.f2[nfft//2-1], FFTop.f2[0],),
                 vmin=0, vmax=np.abs(P).max())
axs[0].set_title(r'$P^+(f, k_x)$')
axs[1].set_title(r'$P^-(f, k_x)$')
axs[1].axis('tight')
axs[1].set_ylim(80, 0);
plt.colorbar(im, ax=axs[1]);

Use function

In [11]:
pup_inv, pdown_inv = WavefieldDecomposition(p, vz, nt, nr, dt, dr, 
                                            rho_sep, vel_sep, nffts=(nfft, nfft), 
                                            critical=critical*100., ntaper=ntaper,
                                            scaling= p.max()/vz.max(), kind='inverse',
                                            dtype='complex128', **dict(damp=1e-10, iter_lim=20, show=2))

print(np.linalg.norm(pdown - pdown_inv))
print(np.linalg.norm(pup - pup_inv))


fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
#axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(vz_obl.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$vzobl$')
#axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(pup_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
#axs[0].plot(r[0], direct, 'r', lw=1)
axs[0].set_title(r'$pup$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(pdown_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$pdown$')
#axs[1].plot(r[0], direct, 'r', lw=1)
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr//2], 'r', lw=2, label=r'$p$')
plt.plot(t, vz_obl[nr//2], '--b', lw=2, label=r'$v_z$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.figure(figsize=(14, 3))
plt.plot(t, pup_inv[nr//2], 'r', lw=2, label=r'$p^-$')
plt.plot(t, pdown_inv[nr//2], '--b', lw=2, label=r'$p^+$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   434562 rows  and   434562 cols
damp = 1.00000000000000e-10   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       20
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   1.361e+01  1.361e+01    1.0e+00  9.8e-02
     1  3.72488e-06   3.221e+00  3.221e+00    2.4e-01  6.4e-01   1.4e+00  1.0e+00
     2 -4.90167e-06   1.643e+00  1.643e+00    1.2e-01  4.2e-01   1.7e+00  2.1e+00
     3 -4.90827e-06   9.905e-01  9.905e-01    7.3e-02  3.2e-01   2.0e+00  3.4e+00
     4 -2.74583e-06   6.425e-01  6.425e-01    4.7e-02  2.7e-01   2.2e+00  4.9e+00
     5 -7.91046e-06   4.318e-01  4.318e-01    3.2e-02  2.4e-01   2.4e+00  6.4e+00
     6 -2.99539e-06   2.958e-01  2.958e-01    2.2e-02  2.2e-01   2.6e+00  8.1e+00
     7 -6.02954e-06   2.087e-01  2.087e-01    1.5e-02  1.9e-01   2.8e+00  9.7e+00
     8 -5.05354e-06   1.549e-01  1.549e-01    1.1e-02  1.6e-01   2.9e+00  1.1e+01
     9 -6.40115e-06   1.227e-01  1.227e-01    9.0e-03  1.4e-01   3.1e+00  1.3e+01
    10 -5.48717e-06   1.034e-01  1.034e-01    7.6e-03  1.1e-01   3.2e+00  1.5e+01
    11 -6.22387e-06   9.155e-02  9.155e-02    6.7e-03  8.8e-02   3.4e+00  1.7e+01
    12 -5.83974e-06   8.401e-02  8.401e-02    6.2e-03  7.3e-02   3.5e+00  1.8e+01
    13 -6.40602e-06   7.881e-02  7.881e-02    5.8e-03  6.2e-02   3.6e+00  2.1e+01
    14 -6.49674e-06   7.484e-02  7.484e-02    5.5e-03  5.5e-02   3.7e+00  2.3e+01
    15 -6.16872e-06   7.153e-02  7.153e-02    5.3e-03  5.0e-02   3.8e+00  2.5e+01
    16 -6.73433e-06   6.866e-02  6.866e-02    5.0e-03  4.7e-02   4.0e+00  2.8e+01
    17 -6.16827e-06   6.617e-02  6.617e-02    4.9e-03  4.3e-02   4.1e+00  3.1e+01
    18 -6.84475e-06   6.406e-02  6.406e-02    4.7e-03  3.9e-02   4.2e+00  3.4e+01
    19 -6.63506e-06   6.232e-02  6.232e-02    4.6e-03  3.5e-02   4.3e+00  3.7e+01
    20 -6.64805e-06   6.087e-02  6.087e-02    4.5e-03  3.1e-02   4.4e+00  4.0e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 6.1e-02   anorm = 4.4e+00   arnorm = 8.4e-03
itn   =      20   r2norm = 6.1e-02   acond = 4.0e+01   xnorm  = 1.1e+01
 
0.877757730395036
0.8777577309633208
Out[11]:
<matplotlib.legend.Legend at 0x1c1ad66ac8>

3. Wavefield separation by inversion with sparsity (model in local linear radon domain)

In [33]:
nwin=31
nwins=8
nover=7
npx=101
pxmax = 5e-4
px = np.linspace(-pxmax, pxmax, npx)
dimsd = p.shape
dims = (nwins*npx, dimsd[1])

# sliding window radon with overlap
Op = Radon2D(t, np.linspace(-dr*nwin//2, dr*nwin//2, nwin), px, centeredh=True,
             kind='linear', engine='numba')
Slidop = Sliding2D(Op, dims, dimsd, nwin, nover, tapertype='cosine', design=True)
dottest(Slidop, np.prod(dimsd), np.prod(dims))

p_pw = Slidop.H * p.ravel()
p_pw = p_pw.reshape(npx*nwins, nt)

fig, axs = plt.subplots(1, 2, figsize=(15, 6))
axs[0].imshow(p.T, cmap='gray',  vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(p_pw.T, cmap='gray',  vmin=-1, vmax=1,
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$p_{pw}$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0);
WARNING: 8 windows required...
WARNING: model wins - start:[  0 101 202 303 404 505 606 707], end:[101 202 303 404 505 606 707 808]
WARNING: data wins - start:[  0  24  48  72  96 120 144 168], end:[ 31  55  79 103 127 151 175 199]
In [ ]:
# subsampling locations
#perc_subsampling = 0.6
#Nsub = int(np.round(par['nx']*perc_subsampling))
#iava = np.sort(np.random.permutation(np.arange(par['nx']))[:Nsub])
iava = np.arange(0, nr, 5)
Nsub= len(iava)

Rop = Restriction(nr*nt, iava, 
                  dims=(nr, nt), 
                  dir=0, dtype='float64')
Rop = BlockDiag([Rop, Rop])


d = np.concatenate((p.flatten(), vz.flatten()))
d_subsampled = Rop*d.flatten()
d_subsampled_adj = Rop.H*d_subsampled.flatten()

d_subsampled = d_subsampled.reshape(2*Nsub, nt)
d_subsampled_adj = d_subsampled_adj.reshape(2*nr, nt)

p_subsampled, vz_subsampled = d_subsampled[:Nsub], d_subsampled[Nsub:]
p_subsampled_adj, vz_subsampled_adj= d_subsampled_adj[:nr], d_subsampled_adj[nr:]


fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(p_subsampled_adj.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'masked $p$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

P = (FFTop*p.flatten()).reshape(nfft, nfft)
PSUB = (FFTop*p_subsampled_adj.flatten()).reshape(nfft, nfft)

PUP_lop = (FFTop*pup_lop.flatten()).reshape(nfft, nfft)
PDOWN_lop = (FFTop*pdown_lop.flatten()).reshape(nfft, nfft)

fig, axs = plt.subplots(1, 2, figsize=(19,5))
im=axs[0].imshow(np.fft.fftshift(np.abs(P[:, :nfft//2-1]),axes=0).T, 
                 cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1],
                           FFTop.f2[nfft//2-1], FFTop.f2[0],),
                 vmin=0, vmax=np.abs(P).max())
axs[0].set_title(r'$P(f, k_x)$')
axs[0].axis('tight')
axs[0].set_ylim(80, 0);
plt.colorbar(im, ax=axs[0])
im=axs[1].imshow(np.fft.fftshift(np.abs(PSUB[:, :nfft//2-1]),axes=0).T,
                 cmap='jet', interpolation='sinc',
                 extent = (FFTop.f1[0], FFTop.f1[nfft//2-1],
                           FFTop.f2[nfft//2-1], FFTop.f2[0],),
                 vmin=0, vmax=np.abs(PSUB).max())
axs[1].set_title(r'$P_{sub}(f, k_x)$')
axs[1].axis('tight')
axs[1].set_ylim(80, 0)
plt.colorbar(im, ax=axs[1]);
In [ ]:
pup_inv, pdown_inv = WavefieldDecomposition(p_subsampled, vz_subsampled, nt, nr, dt, dr,
                                            rho_sep, vel_sep, 
                                            nffts=(nfft, nfft), scaling= p.max()/vz.max(),
                                            sptransf=Slidop, restriction=Rop,
                                            dtype='complex128', **dict(damp=1e-10, iter_lim=50, show=2))

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(vz_obl.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[1].set_title(r'$vzobl$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)


fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(pup_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[0].set_title(r'$pup$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(pdown_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[1].set_title(r'$pdown$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr//2], 'r', lw=2, label=r'$p$')
plt.plot(t, vz_obl[nr//2], '--b', lw=2, label=r'$v_z$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.figure(figsize=(14, 3))
plt.plot(t, pup_inv[nr//2], 'r', lw=2, label=r'$p^-$')
plt.plot(t, pdown_inv[nr//2], '--b', lw=2, label=r'$p^+$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend();
 
LSQR            Least-squares solution of  Ax = b
The matrix A has    88642 rows  and 1.7469e+06 cols
damp = 1.00000000000000e-10   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       50
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   6.103e+00  6.103e+00    1.0e+00  1.8e+00
     1 -8.01936e-09   2.186e+00  2.186e+00    3.6e-01  8.4e-01   1.2e+01  1.0e+00
     2  6.04077e-07   9.743e-01  9.743e-01    1.6e-01  7.0e-01   1.7e+01  2.1e+00
     3 -8.50818e-07   7.420e-01  7.420e-01    1.2e-01  4.9e-01   2.6e+01  3.6e+00
     4 -9.97581e-08   5.325e-01  5.325e-01    8.7e-02  2.8e-01   3.5e+01  5.5e+00
     5  4.43181e-07   4.654e-01  4.654e-01    7.6e-02  1.8e-01   4.1e+01  7.0e+00
     6 -3.12334e-07   4.202e-01  4.202e-01    6.9e-02  1.2e-01   4.6e+01  8.7e+00
     7 -3.51240e-07   4.039e-01  4.039e-01    6.6e-02  6.6e-02   5.2e+01  1.0e+01
     8  2.39867e-07   3.961e-01  3.961e-01    6.5e-02  5.0e-02   5.6e+01  1.2e+01
     9  3.56111e-07   3.916e-01  3.916e-01    6.4e-02  3.4e-02   6.1e+01  1.4e+01
    10 -1.81410e-07   3.892e-01  3.892e-01    6.4e-02  3.1e-02   6.5e+01  1.5e+01
In [ ]:
pup_inv, pdown_inv = WavefieldDecomposition(p_subsampled, vz_subsampled, nt, nr, dt, dr,
                                            rho_sep, vel_sep, 
                                            nffts=(nfft, nfft), scaling= p.max()/vz.max(),
                                            critical=critical*100, ntaper=ntaper,
                                            sptransf=Slidop, restriction=Rop, solver=FISTA,
                                            dtype='complex128', **dict(niter=40, eps=8e-2, show=True))

fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(p_subsampled_adj.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow((p.max()/vz.max())*vz_subsampled_adj.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[1].set_title(r'$vzobl$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)


fig, axs = plt.subplots(1, 2, figsize=(9,6))
axs[0].imshow(pup_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[0].set_title(r'$pup$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(pdown_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (x.min(),x.max(),t.max(),t.min()))
axs[1].set_title(r'$pdown$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr//2], 'r', lw=2, label=r'$p$')
plt.plot(t, vz_obl[nr//2], '--b', lw=2, label=r'$v_z$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.figure(figsize=(14, 3))
plt.plot(t, pup_inv[nr//2], 'r', lw=2, label=r'$p^-$')
plt.plot(t, pdown_inv[nr//2], '--b', lw=2, label=r'$p^+$')
plt.xlim(0,t[200])
plt.ylim(-0.2, 0.2)
plt.legend();
In [ ]:
p_inv_pw = Slidop.H * (pup_inv.ravel() + pdown_inv.ravel())
p_inv_pw = p_inv_pw.reshape(npx*nwins, nt)

fig, axs = plt.subplots(1, 4, figsize=(15, 6))
axs[0].imshow(p.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[0].set_title(r'$p$')
axs[0].axis('tight')
axs[0].set_ylim(1, 0)
axs[1].imshow(p_pw.T, cmap='gray',  vmin=-1, vmax=1,
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[1].set_title(r'$p_{pw} rec$')
axs[1].axis('tight')
axs[1].set_ylim(1, 0);
axs[2].imshow(pup_inv.T + pdown_inv.T, cmap='gray', vmin=-0.1*np.abs(p).max(), vmax=0.1*np.abs(p).max(),
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[2].set_title(r'$p^{inv}_{pw}$')
axs[2].axis('tight')
axs[2].set_ylim(1, 0);
axs[3].imshow(p_inv_pw.T, cmap='gray',  vmin=-1, vmax=1,
              extent = (r[0, 0], r[0,-1], t[-1], t[0]))
axs[3].set_title(r'$p^{inv}_{pw}$')
axs[3].axis('tight')
axs[3].set_ylim(1, 0);